Univariate Arima Models

Predicting COVID-19 Deaths (STAT 390)

Author

Ryan Nguyen

1. Load Libraries, Seed, & Data

Code
# Load necessary libraries
library(forecast)
library(tseries)
library(tidyverse)
library(tidymodels)
library(kableExtra)
library(modeltime)
library(timetk)
library(cowplot)

# setting the seed
set.seed(123)

# Loading Data  ----------------------------------------------------------------
daily_deaths <- read_csv("data/new_daily_deaths.csv") %>% 
  select(daily_deaths)

east <- read_csv("data/east_daily.csv")

midwest <- read_csv("data/midwest_daily.csv")

south <- read_csv("data/south_daily.csv")

west <- read_csv("data/west_daily.csv")

2. Data Splitting

Code
# Data Splitting ---------------------------------------------------------------
east_splits <- time_series_split(east, initial = "1022 days", assess = "114 days")
east_train <- training(east_splits)
east_test <- testing(east_splits)

midwest_splits <- time_series_split(midwest, initial = "1022 days", assess = "114 days")
midwest_train <- training(midwest_splits)
midwest_test <- testing(midwest_splits)

south_splits <- time_series_split(south, initial = "1022 days", assess = "114 days")
south_train <- training(south_splits)
south_test <- testing(south_splits)

west_splits <- time_series_split(west, initial = "1022 days", assess = "114 days")
west_train <- training(west_splits)
west_test <- testing(west_splits)

## Convert Into Time Series ----------------------------------------------------
east_train_ts <- ts(east_train %>%
                      select(daily_deaths))
east_test_ts <- ts(east_test %>%
                     select(daily_deaths))
east_ts <- ts(east %>% 
                select(daily_deaths))

midwest_train_ts <- ts(midwest_train %>%
                      select(daily_deaths))
midwest_test_ts <- ts(midwest_test %>%
                     select(daily_deaths))
midwest_ts <- ts(midwest %>% 
                select(daily_deaths))

south_train_ts <- ts(south_train %>%
                    select(daily_deaths))
south_test_ts <- ts(south_test %>%
                   select(daily_deaths))
south_ts <- ts(south %>%
                 select(daily_deaths))

west_train_ts <- ts(west_train %>%
                   select(daily_deaths))
west_test_ts <- ts(west_test %>%
                  select(daily_deaths))
west_ts <- ts(west %>%
                select(daily_deaths))

## Plotting Training and Testing -----------------------------------------------
par(mfrow=c(2,2))
ggplot() +
  geom_line(data = east_test, aes(x = date, y = daily_deaths, color = "turquoise")) +
  geom_point(data = east_test, aes(x = date, y = daily_deaths, color = "turquoise"), size = .5) +
  geom_line(data = east_train, aes(x = date, y = daily_deaths, color = "orange")) +
  geom_point(data = east_train, aes(x = date, y = daily_deaths, color = "orange"), size = .5) +
  scale_color_manual(name="Data Set", values = c("turquoise", "orange"), labels = c("Training", "Testing")) +
  labs(y = "Daily Deaths",
       x = "Date",
       title = "East Region Training vs. Testing Data")

Code
ggplot() +
  geom_line(data = midwest_test, aes(x = date, y = daily_deaths, color = "turquoise")) +
  geom_point(data = midwest_test, aes(x = date, y = daily_deaths, color = "turquoise"), size = .5) +
  geom_line(data = midwest_train, aes(x = date, y = daily_deaths, color = "orange")) +
  geom_point(data = midwest_train, aes(x = date, y = daily_deaths, color = "orange"), size = .5) +
  scale_color_manual(name="Data Set", values = c("turquoise", "orange"), labels = c("Training", "Testing")) +
  labs(y = "Daily Deaths",
       x = "Date",
       title = "Midwest Region Training vs. Testing Data")

Code
ggplot() +
  geom_line(data = south_test, aes(x = date, y = daily_deaths, color = "turquoise")) +
  geom_point(data = south_test, aes(x = date, y = daily_deaths, color = "turquoise"), size = .5) +
  geom_line(data = south_train, aes(x = date, y = daily_deaths, color = "orange")) +
  geom_point(data = south_train, aes(x = date, y = daily_deaths, color = "orange"), size = .5) +
  scale_color_manual(name="Data Set", values = c("turquoise", "orange"), labels = c("Training", "Testing")) +
  labs(y = "Daily Deaths",
       x = "Date",
       title = "South Region Training vs. Testing Data")

Code
ggplot() +
  geom_line(data = west_test, aes(x = date, y = daily_deaths, color = "turquoise")) +
  geom_point(data = west_test, aes(x = date, y = daily_deaths, color = "turquoise"), size = .5) +
  geom_line(data = west_train, aes(x = date, y = daily_deaths, color = "orange")) +
  geom_point(data = west_train, aes(x = date, y = daily_deaths, color = "orange"), size = .5) +
  scale_color_manual(name="Data Set", values = c("turquoise", "orange"), labels = c("Training", "Testing")) +
  labs(y = "Daily Deaths",
       x = "Date",
       title = "West Region Training vs. Testing Data")

3. Testing for Stationarity

Code
### Using Augmented Dickey-Fuller Test -----------------------------------------
adf_east <- adf.test(east_ts)
print(adf_east) # p-value is 0.01 < 0.05, implying it is stationary

    Augmented Dickey-Fuller Test

data:  east_ts
Dickey-Fuller = -4.5566, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
Code
adf_midwest <- adf.test(midwest_ts)
print(adf_midwest) # p-value is 0.4919 > 0.05, implying it is NOT stationary

    Augmented Dickey-Fuller Test

data:  midwest_ts
Dickey-Fuller = -2.5984, Lag order = 10, p-value = 0.325
alternative hypothesis: stationary
Code
midwest_ts <- diff(midwest_ts)
adf_midwest <- adf.test(midwest_ts)
print(adf_midwest) # p-value is 0.01 < 0.05, implying it is stationary with 1 diff

    Augmented Dickey-Fuller Test

data:  midwest_ts
Dickey-Fuller = -13.43, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
Code
adf_south <- adf.test(south_ts)
print(adf_south) # p-value is 0.3018 > 0.05, implying it is NOT stationary

    Augmented Dickey-Fuller Test

data:  south_ts
Dickey-Fuller = -2.6533, Lag order = 10, p-value = 0.3018
alternative hypothesis: stationary
Code
south_ts <- diff(south_ts)
adf_south <- adf.test(south_ts)
print(adf_south) # p-value is 0.01 < 0.05, implying it is stationary with 1 diff

    Augmented Dickey-Fuller Test

data:  south_ts
Dickey-Fuller = -11.36, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
Code
adf_west <- adf.test(west_ts)
print(adf_west) # p-value is 0.4098 > 0.05, implying it is NOT stationary

    Augmented Dickey-Fuller Test

data:  west_ts
Dickey-Fuller = -2.3981, Lag order = 10, p-value = 0.4098
alternative hypothesis: stationary
Code
west_ts <- diff(west_ts)
adf_west <- adf.test(west_ts)
print(adf_west)  # p-value is 0.01 < 0.05, implying it is stationary with 1 diff

    Augmented Dickey-Fuller Test

data:  west_ts
Dickey-Fuller = -12.247, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
Code
### Making Testing/Training Stationary -----------------------------------------
midwest_train_ts <- diff(midwest_train_ts)
midwest_test_ts <- diff(midwest_test_ts)

south_train_ts <- diff(south_train_ts)
south_test_ts <- diff(south_test_ts)

west_train_ts <- diff(west_train_ts)
west_test_ts <- diff(west_test_ts)

4. ACF & PACF Plots

Code
## ACF & PACF Plots ------------------------------------------------------------
par(mfrow=c(3, 1))
acf(east_ts, main = "East")
pacf(east_ts, main = "East")
plot(east_ts, main = "Daily Deaths: East")

Code
acf(midwest_ts, main = "Midwest")
pacf(midwest_ts, main = "Midwest")
plot(midwest_ts, main = "Daily Deaths: Midwest")

Code
acf(south_ts, main = "South")
pacf(south_ts, main = "South")
plot(south_ts, main = "Daily Deaths: South")

Code
acf(west_ts, main = "West")
pacf(west_ts, main = "West")
plot(south_ts, main = "Daily Deaths: South")

Code
plot(west_ts, main = "Daily Deaths: West")

5. Time Series Decomposition

Code
# total daily death decomposition
daily_ts <- ts(daily_deaths, frequency = 365)

decomposed_ts <- decompose(daily_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(daily_ts, main = "Total Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_ts$trend, main = "Trend Component", ylab = "Value", col = "red")

Code
plot(decomposed_ts$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_ts$random, main = "Residual Component", ylab = "Value", col = "purple")

Code
# East Decomposition -----------------------------------------------------------
east_daily <- east %>% 
  select(daily_deaths)

east_ts <- ts(east_daily, frequency = 365)

decomposed_east <- decompose(east_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(east_ts, main = "East Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_east$trend, main = "Trend Component", ylab = "Value", col = "red")

Code
plot(decomposed_east$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_east$random, main = "Residual Component", ylab = "Value", col = "purple")

Code
## Midwest Decomposition ----------------------------------------------------------
midwest_daily <- midwest %>% 
  select(daily_deaths)

midwest_ts <- ts(midwest_daily, frequency = 365)

decomposed_midwest <- decompose(midwest_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(midwest_ts, main = "Midwest Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_midwest$trend, main = "Trend Component", ylab = "Value", col = "red")

Code
plot(decomposed_midwest$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_midwest$random, main = "Residual Component", ylab = "Value", col = "purple")

Code
## South Decomposition ----------------------------------------------------------
south_daily <- south %>% 
  select(daily_deaths)

south_ts <- ts(south_daily, frequency = 365)

decomposed_south <- decompose(south_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(south_ts, main = "South Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_south$trend, main = "Trend Component", ylab = "Value", col = "red")

Code
plot(decomposed_south$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_south$random, main = "Residual Component", ylab = "Value", col = "purple")

Code
## West Decomposition ----------------------------------------------------------
west_daily <- west %>% 
  select(daily_deaths)

west_ts <- ts(west_daily, frequency = 365)

decomposed_west <- decompose(west_ts)

# Plot the original time series and the decomposed components
par(mfrow=c(2,1))
plot(west_ts, main = "West Original Time Series", ylab = "Value", col = "blue")
plot(decomposed_west$trend, main = "Trend Component", ylab = "Value", col = "red")

Code
plot(decomposed_west$seasonal, main = "Seasonal Component", ylab = "Value", col = "green")
plot(decomposed_west$random, main = "Residual Component", ylab = "Value", col = "purple")

6. ARIMA Model

6.1 Initial Model

Code
# Building ARIMA Model ---------------------------------------------------------
east_arima <- arima(east_train_ts, order=c(1,0,1))
midwest_arima <- arima(midwest_train_ts, order=c(1,1,1))
south_arima <- arima(south_train_ts, order=c(1,1,1))
west_arima <- arima(west_train_ts, order=c(1,1,1))

## Summary & Forecast ---------------------------------------------------------
summary(east_arima)

Call:
arima(x = east_train_ts, order = c(1, 0, 1))

Coefficients:
         ar1      ma1  intercept
      0.9176  -0.0388   197.3738
s.e.  0.0147   0.0493    39.7998

sigma^2 estimated as 12153:  log likelihood = -6257.17,  aic = 12522.33

Training set error measures:
                    ME     RMSE      MAE  MPE MAPE      MASE        ACF1
Training set 0.1151898 110.2407 64.66739 -Inf  Inf 0.9918341 0.008515868
Code
summary(midwest_arima)

Call:
arima(x = midwest_train_ts, order = c(1, 1, 1))

Coefficients:
          ar1      ma1
      -0.3477  -1.0000
s.e.   0.0293   0.0026

sigma^2 estimated as 50173:  log likelihood = -6970.99,  aic = 13947.98

Training set error measures:
                     ME     RMSE      MAE MPE MAPE      MASE        ACF1
Training set -0.8748943 223.8834 138.9352 NaN  Inf 0.6215031 -0.06581815
Code
summary(south_arima)

Call:
arima(x = south_train_ts, order = c(1, 1, 1))

Coefficients:
          ar1      ma1
      -0.1727  -1.0000
s.e.   0.0308   0.0026

sigma^2 estimated as 106737:  log likelihood = -7355.8,  aic = 14717.6

Training set error measures:
                    ME    RMSE      MAE MPE MAPE      MASE        ACF1
Training set -1.420592 326.547 203.7282 NaN  Inf 0.6344327 -0.04113314
Code
summary(west_arima)

Call:
arima(x = west_train_ts, order = c(1, 1, 1))

Coefficients:
          ar1      ma1
      -0.1698  -1.0000
s.e.   0.0309   0.0026

sigma^2 estimated as 23955:  log likelihood = -6593.77,  aic = 13193.54

Training set error measures:
                     ME     RMSE      MAE MPE MAPE      MASE        ACF1
Training set -0.6839532 154.6997 101.0808 NaN  Inf 0.6424269 -0.02682684
Code
east_values <- forecast(east_arima, h = length(east_test_ts))
midwest_values <- forecast(midwest_arima, h = length(midwest_test_ts))
south_values <- forecast(south_arima, h = length(south_test_ts))
west_values <- forecast(west_arima, h = length(west_test_ts))

6.2 Accuracy

Code
## Evaluating Accuracy ---------------------------------------------------------
east_initial_arima_acc <- forecast::accuracy(east_values)

east_initial_arima_acc %>% 
  kbl(caption = "<center>Initial East Accuracy<center>") %>% 
  kable_styling()
Initial East Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1151898 110.2407 64.66739 -Inf Inf 0.9918341 0.0085159
Code
### East Region
east_test_plot <- as.data.frame(east_test_ts) %>% 
  mutate(num = seq(1023, (1023+113)))

east_initial_arima_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Initial ARIMA",
       y = "Daily Deaths")

east_initial_arima_plot

Code
### Midwest Region

midwest_initial_arima_acc <- forecast::accuracy(midwest_values)

midwest_initial_arima_acc %>% 
  kbl(caption = "<center>Initial Midwest Accuracy<center>") %>% 
  kable_styling()
Initial Midwest Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.8748943 223.8834 138.9352 NaN Inf 0.6215031 -0.0658182
Code
midwest_test_plot <- as.data.frame(midwest_test_ts) %>% 
  mutate(num = seq(1023, (1023+112)))

midwest_initial_arima_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Initial ARIMA",
       y = "Daily Deaths")

midwest_initial_arima_plot

Code
### South Region

south_initial_arima_acc <- forecast::accuracy(south_values)

south_initial_arima_acc %>% 
  kbl(caption = "<center>Initial South Accuracy<center>") %>% 
  kable_styling()
Initial South Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set -1.420592 326.547 203.7282 NaN Inf 0.6344327 -0.0411331
Code
south_test_plot <- as.data.frame(south_test_ts) %>% 
  mutate(num = seq(1023, (1023+112)))

south_initial_arima_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Initial ARIMA",
       y = "Daily Deaths")

south_initial_arima_plot 

Code
### West Region

west_initial_arima_acc <- forecast::accuracy(west_values)

west_initial_arima_acc %>% 
  kbl(caption = "<center>Initial West Accuracy<center>") %>% 
  kable_styling()
Initial West Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.6839532 154.6997 101.0808 NaN Inf 0.6424269 -0.0268268
Code
west_test_plot <- as.data.frame(west_test_ts) %>% 
  mutate(num = seq(1023, (1023+112)))

west_initial_arima_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Initial ARIMA",
       y = "Daily Deaths") 

west_initial_arima_plot

6.4 Rebuilding ARIMA Model

Code
# Rebuilding ARIMA Models ------------------------------------------------------
east_arima <- arima(east_train_ts, order=c(5,0,5))
midwest_arima <- arima(midwest_train_ts, order=c(5,1,5))
south_arima <- arima(south_train_ts, order=c(5,1,5))
west_arima <- arima(west_train_ts, order=c(5,1,5), method = "ML")


east_values <- forecast(east_arima, h = length(east_test_ts))
midwest_values <- forecast(midwest_arima, h = length(midwest_test_ts))
south_values <- forecast(south_arima, h = length(south_test_ts))
west_values <- forecast(west_arima, h = length(west_test_ts))

## Evaluating Accuracy & Plots -------------------------------------------------
east_tuned_arima_acc <- forecast::accuracy(east_values)

east_tuned_arima_acc %>% 
  kbl(caption = "<center>Tuned East Accuracy<center>") %>% 
  kable_styling()
Tuned East Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1095068 80.35755 47.18823 NaN Inf 0.723748 -0.0254917
Code
### East Region
east_tuned_arima_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Tuned ARIMA",
       y = "Daily Deaths")

east_tuned_arima_plot

Code
### Midwest Region

midwest_tuned_arima_acc <- forecast::accuracy(midwest_values)

midwest_tuned_arima_acc %>% 
  kbl(caption = "<center>Tuned Midwest Accuracy<center>") %>% 
  kable_styling()
Tuned Midwest Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set -3.635509 152.8728 95.00769 NaN Inf 0.4250009 -0.0128029
Code
midwest_tuned_arima_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Tuned ARIMA",
       y = "Daily Deaths")

midwest_tuned_arima_plot

Code
### South Region

south_tuned_arima_acc <- forecast::accuracy(south_values)

south_tuned_arima_acc %>% 
  kbl(caption = "<center>Tuned South Accuracy<center>") %>% 
  kable_styling()
Tuned South Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set -4.199063 200.7617 120.224 NaN Inf 0.3743912 -0.0666043
Code
south_tuned_arima_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Tuned ARIMA",
       y = "Daily Deaths")

south_tuned_arima_plot 

Code
### West Region
west_tuned_arima_acc <- forecast::accuracy(west_values)

west_tuned_arima_acc %>% 
  kbl(caption = "<center>Tuned West Accuracy<center>") %>% 
  kable_styling()
Tuned West Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.0352798 99.27405 63.93076 NaN Inf 0.4063168 -0.0093102
Code
west_tuned_arima_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Tuned ARIMA",
       y = "Daily Deaths") 

west_tuned_arima_plot

Code
# saving results
save(east_tuned_arima_acc, east_initial_arima_acc, east_initial_arima_plot, east_tuned_arima_plot,
     midwest_tuned_arima_acc, midwest_initial_arima_acc, midwest_initial_arima_plot, midwest_tuned_arima_plot,
     south_tuned_arima_acc, south_initial_arima_acc, south_initial_arima_plot, south_tuned_arima_plot,
     west_tuned_arima_acc, west_initial_arima_acc, west_initial_arima_plot, west_tuned_arima_plot,
     file = "results/uni_arima.rda"
     )

7. SARIMA Model

7.1 Initial SARIMA Model

Code
# Building SARIMA model  ------------------------------------------------------
east_sarima <- arima(east_train_ts, seasonal = list(order = c(1,1,1)))
midwest_sarima <- arima(midwest_train_ts, seasonal = list(order = c(1,1,1)))
south_sarima <- arima(south_train_ts, seasonal = list(order = c(1,1,1)))
west_sarima <- arima(west_train_ts, seasonal = list(order = c(1,1,1)))

# Summary of the SARIMA model
summary(east_sarima)

Call:
arima(x = east_train_ts, seasonal = list(order = c(1, 1, 1)))

Coefficients:
        sar1     sma1
      0.4706  -0.8016
s.e.  0.0382   0.0211

sigma^2 estimated as 11156:  log likelihood = -6206.64,  aic = 12419.27

Training set error measures:
                    ME     RMSE     MAE  MPE MAPE      MASE      ACF1
Training set 0.1606608 105.5715 61.2767 -Inf  Inf 0.9398295 0.1140208
Code
summary(midwest_sarima)

Call:
arima(x = midwest_train_ts, seasonal = list(order = c(1, 1, 1)))

Coefficients:
         sar1     sma1
      -0.3477  -1.0000
s.e.   0.0293   0.0026

sigma^2 estimated as 50173:  log likelihood = -6970.99,  aic = 13947.98

Training set error measures:
                     ME     RMSE      MAE MPE MAPE      MASE        ACF1
Training set -0.8748943 223.8834 138.9352 NaN  Inf 0.6215031 -0.06581815
Code
summary(south_sarima)

Call:
arima(x = south_train_ts, seasonal = list(order = c(1, 1, 1)))

Coefficients:
         sar1     sma1
      -0.1727  -1.0000
s.e.   0.0308   0.0026

sigma^2 estimated as 106737:  log likelihood = -7355.8,  aic = 14717.6

Training set error measures:
                    ME    RMSE      MAE MPE MAPE      MASE        ACF1
Training set -1.420592 326.547 203.7282 NaN  Inf 0.6344327 -0.04113314
Code
summary(west_sarima)

Call:
arima(x = west_train_ts, seasonal = list(order = c(1, 1, 1)))

Coefficients:
         sar1     sma1
      -0.1698  -1.0000
s.e.   0.0309   0.0026

sigma^2 estimated as 23955:  log likelihood = -6593.77,  aic = 13193.54

Training set error measures:
                     ME     RMSE      MAE MPE MAPE      MASE        ACF1
Training set -0.6839532 154.6997 101.0808 NaN  Inf 0.6424269 -0.02682684
Code
## Forecasting with the SARIMA model -------------------------------------------
east_values <- forecast(east_sarima, h = length(east_test_ts))
midwest_values <- forecast(midwest_sarima, h = length(midwest_test_ts))
south_values <- forecast(south_sarima, h = length(south_test_ts))
west_values <- forecast(west_sarima, h = length(west_test_ts))

7.2 Accuracy

Code
## Plotting & Accuracy

## East Region
east_initial_sarima_acc <- forecast::accuracy(east_values)

east_initial_sarima_acc %>% 
  kbl(caption = "<center>Initial SARIMA East Accuracy<center>") %>% 
  kable_styling()
Initial SARIMA East Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1606608 105.5715 61.2767 -Inf Inf 0.9398295 0.1140208
Code
east_initial_sarima_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Initial SARIMA",
       y = "Daily Deaths")

east_initial_sarima_plot

Code
## Midwest Region
midwest_initial_sarima_acc <- forecast::accuracy(midwest_values)

midwest_initial_sarima_acc %>% 
  kbl(caption = "<center>Initial SARIMA Midwest Accuracy<center>") %>% 
  kable_styling()
Initial SARIMA Midwest Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.8748943 223.8834 138.9352 NaN Inf 0.6215031 -0.0658182
Code
midwest_initial_sarima_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Initial SARIMA",
       y = "Daily Deaths")

midwest_initial_sarima_plot

Code
## South Region
south_initial_sarima_acc <- forecast::accuracy(south_values)

south_initial_sarima_acc %>% 
  kbl(caption = "<center>Initial SARIMA South Accuracy<center>") %>% 
  kable_styling()
Initial SARIMA South Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set -1.420592 326.547 203.7282 NaN Inf 0.6344327 -0.0411331
Code
south_initial_sarima_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Initial SARIMA",
       y = "Daily Deaths")

south_initial_sarima_plot

Code
## West Region
west_initial_sarima_acc <- forecast::accuracy(west_values)

west_initial_sarima_acc %>% 
  kbl(caption = "<center>Initial SARIMA West Accuracy<center>") %>% 
  kable_styling()
Initial SARIMA West Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.6839532 154.6997 101.0808 NaN Inf 0.6424269 -0.0268268
Code
west_initial_sarima_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Initial SARIMA",
       y = "Daily Deaths")

west_initial_sarima_plot

7.3 Grid Search

Code
## Define the parameter grids --------------------------------------------------
p_grid <- 0:5  # AR parameter
d_grid <- 0:0  # Differencing parameter
q_grid <- 0:5  # MA parameter
P_grid <- 0:1  # Seasonal AR parameter
D_grid <- 0:1  # Seasonal differencing parameter
Q_grid <- 0:1  # Seasonal MA parameter
s <- 12        # Seasonal period



## East Grid Search ------------------------------------------------------------

# Create East Region
east_train <- east_train %>% 
  select(daily_deaths)

# Create an empty data frame to store the results
east_results <- data.frame(order = character(),
                      seasonal = character(),
                      AIC = numeric(),
                      stringsAsFactors = FALSE)

# Perform grid search
for (p in p_grid) {
  for (d in d_grid) {
    for (q in q_grid) {
      for (P in P_grid) {
        for (D in D_grid) {
          for (Q in Q_grid) {
            order <- c(p, d, q)
            seasonal <- c(P, D, Q, s)
            model <- tryCatch(
              {
                fit <- arima(east_train, order = order, seasonal = seasonal)
                AIC_val <- AIC(fit)
                # Append the results to the data frame
                east_results <- rbind(east_results, data.frame(order = paste(order, collapse = ", "),
                                                     seasonal = paste(seasonal, collapse = ", "),
                                                     AIC = AIC_val))
              },
              error = function(e) NULL
            )
          }
        }
      }
    }
  }
}

d_grid <- 0:1  # Differencing parameter

## Midwest Grid Search ------------------------------------------------------------

# Create Midwest Region
midwest_train <- midwest_train %>% 
  select(daily_deaths)

# Create an empty data frame to store the results
midwest_results <- data.frame(order = character(),
                           seasonal = character(),
                           AIC = numeric(),
                           stringsAsFactors = FALSE)

# Perform grid search
for (p in p_grid) {
  for (d in d_grid) {
    for (q in q_grid) {
      for (P in P_grid) {
        for (D in D_grid) {
          for (Q in Q_grid) {
            order <- c(p, d, q)
            seasonal <- c(P, D, Q, s)
            model <- tryCatch(
              {
                fit <- arima(midwest_train, order = order, seasonal = seasonal)
                AIC_val <- AIC(fit)
                # Append the results to the data frame
                midwest_results <- rbind(midwest_results, data.frame(order = paste(order, collapse = ", "),
                                                               seasonal = paste(seasonal, collapse = ", "),
                                                               AIC = AIC_val))
              },
              error = function(e) NULL
            )
          }
        }
      }
    }
  }
}

## South Grid Search ------------------------------------------------------------

# Create South Region
south_train <- south_train %>% 
  select(daily_deaths)

# Create an empty data frame to store the results
south_results <- data.frame(order = character(),
                           seasonal = character(),
                           AIC = numeric(),
                           stringsAsFactors = FALSE)

# Perform grid search
for (p in p_grid) {
  for (d in d_grid) {
    for (q in q_grid) {
      for (P in P_grid) {
        for (D in D_grid) {
          for (Q in Q_grid) {
            order <- c(p, d, q)
            seasonal <- c(P, D, Q, s)
            model <- tryCatch(
              {
                fit <- arima(south_train, order = order, seasonal = seasonal)
                AIC_val <- AIC(fit)
                # Append the results to the data frame
                south_results <- rbind(south_results, data.frame(order = paste(order, collapse = ", "),
                                                               seasonal = paste(seasonal, collapse = ", "),
                                                               AIC = AIC_val))
              },
              error = function(e) NULL
            )
          }
        }
      }
    }
  }
}


## West Grid Search ------------------------------------------------------------

# Create West Region
west_train <- west_train %>% 
  select(daily_deaths)

# Create an empty data frame to store the results
west_results <- data.frame(order = character(),
                           seasonal = character(),
                           AIC = numeric(),
                           stringsAsFactors = FALSE)

# Perform grid search
for (p in p_grid) {
  for (d in d_grid) {
    for (q in q_grid) {
      for (P in P_grid) {
        for (D in D_grid) {
          for (Q in Q_grid) {
            order <- c(p, d, q)
            seasonal <- c(P, D, Q, s)
            model <- tryCatch(
              {
                fit <- arima(west_train, order = order, seasonal = seasonal)
                AIC_val <- AIC(fit)
                # Append the results to the data frame
                west_results <- rbind(west_results, data.frame(order = paste(order, collapse = ", "),
                                                               seasonal = paste(seasonal, collapse = ", "),
                                                               AIC = AIC_val))
              },
              error = function(e) NULL
            )
          }
        }
      }
    }
  }
}
Code
save(east_results, midwest_results, south_results, west_results, file = "results/sarima_grid.rda")
Code
load("results/sarima_grid.rda")

Best East SARIMA

Code
# Find the best model (minimum AIC)
best_east_model <- east_results[which.min(east_results$AIC), ]
print(best_east_model)
      order    seasonal      AIC
213 5, 0, 4 0, 1, 1, 12 11865.53

Best Midwest SARIMA

Code
# Midwest Results
# Find the best model (minimum AIC)
best_midwest_model <- midwest_results[which.min(midwest_results$AIC), ]
print(best_midwest_model)
      order    seasonal      AIC
385 5, 1, 5 0, 0, 1, 12 12819.82

Best South SARIMA

Code
# South Results
# Find the best model (minimum AIC)
best_south_model <- south_results[which.min(south_results$AIC), ]
print(best_south_model)
      order    seasonal      AIC
402 5, 1, 5 0, 0, 0, 12 13698.93

Best West SARIMA

Code
# Find the best model (minimum AIC)
best_west_model <- west_results[which.min(west_results$AIC), ]
print(best_west_model)
      order    seasonal      AIC
385 5, 0, 5 1, 1, 0, 12 12141.63

7.4 Rebuilding SARIMA

Code
# Rebuilding SARIMA Models -----------------------------------------------------
east_sarima <- arima(east_train_ts, order = c(5, 0, 4), 
                     seasonal = list(order = c(0, 1, 1), period = 12))
midwest_sarima <- arima(midwest_train_ts, order = c(5, 1, 5),
                        seasonal = list(order = c(0, 0, 1), period = 12))
south_sarima <- arima(south_train_ts, order = c(5, 1, 5),
                           seasonal = list(order = c(0, 0, 0), period = 12))
west_sarima <- arima(west_train_ts, order = c(5, 0, 5),
                     seasonal = list(order = c(1, 1, 0), period = 12))


## Forecasting with the NEW SARIMA model ---------------------------------------
east_values <- forecast(east_sarima, h = length(east_test_ts))
midwest_values <- forecast(midwest_sarima, h = length(midwest_test_ts))
south_values <- forecast(south_sarima, h = length(south_test_ts))
west_values <- forecast(west_sarima, h = length(west_test_ts))
Code
## Plotting & Accuracy

## East Region
east_tuned_sarima_acc <- forecast::accuracy(east_values)

east_tuned_sarima_acc %>% 
  kbl(caption = "<center>Tuned SARIMA East Accuracy<center>") %>% 
  kable_styling()
Tuned SARIMA East Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1031443 87.81972 53.94287 NaN Inf 0.8273471 -0.0998164
Code
east_tuned_sarima_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Tuned SARIMA",
       y = "Daily Deaths")

east_tuned_sarima_plot

Code
## Midwest Region
midwest_tuned_sarima_acc <- forecast::accuracy(midwest_values)

midwest_tuned_sarima_acc %>% 
  kbl(caption = "<center>Tuned SARIMA Midwest Accuracy<center>") %>% 
  kable_styling()
Tuned SARIMA Midwest Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.1147262 143.4927 91.79687 NaN Inf 0.4106378 -0.012698
Code
midwest_tuned_sarima_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Tuned SARIMA",
       y = "Daily Deaths")

midwest_tuned_sarima_plot

Code
## South Region
south_tuned_sarima_acc <- forecast::accuracy(south_values)

south_tuned_sarima_acc %>% 
  kbl(caption = "<center>Tuned SARIMA South Accuracy<center>") %>% 
  kable_styling()
Tuned SARIMA South Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set -4.199063 200.7617 120.224 NaN Inf 0.3743912 -0.0666043
Code
south_tuned_sarima_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Tuned SARIMA",
       y = "Daily Deaths")

south_tuned_sarima_plot

Code
## West Region
west_tuned_sarima_acc <- forecast::accuracy(west_values)

west_tuned_sarima_acc %>% 
  kbl(caption = "<center>Tuned SARIMA West Accuracy<center>") %>% 
  kable_styling()
Tuned SARIMA West Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.0032713 108.9033 71.27948 NaN Inf 0.4530221 0.0086743
Code
west_tuned_sarima_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Tuned SARIMA",
       y = "Daily Deaths")

west_tuned_sarima_plot

Code
# saving results
save(east_tuned_sarima_acc, east_initial_sarima_acc, east_initial_sarima_plot, east_tuned_sarima_plot,
     midwest_tuned_sarima_acc, midwest_initial_sarima_acc, midwest_initial_sarima_plot, midwest_tuned_sarima_plot,
     south_tuned_sarima_acc, south_initial_sarima_acc, south_initial_sarima_plot, south_tuned_sarima_plot,
     west_tuned_sarima_acc, west_initial_sarima_acc, west_initial_sarima_plot, west_tuned_sarima_plot,
     file = "results/uni_sarima.rda"
     )

8. Auto ARIMA Model

8.1 Building AutoARIMA

Code
# Building AutoARIMA Model -----------------------------------------------------
east_arima <- auto.arima(east_train_ts, seasonal = TRUE)
midwest_arima <- auto.arima(midwest_train_ts, seasonal = TRUE)
south_arima <- auto.arima(south_train_ts, seasonal = TRUE)
west_arima <- auto.arima(west_train_ts, seasonal = TRUE)

# Summary of the ARIMA model
summary(east_arima)
Series: east_train_ts 
ARIMA(5,1,2) 

Coefficients:
         ar1      ar2      ar3      ar4      ar5      ma1     ma2
      0.5142  -0.7304  -0.1394  -0.3479  -0.2207  -1.0295  0.7353
s.e.  0.0480   0.0344   0.0426   0.0329   0.0418   0.0410  0.0265

sigma^2 = 6861:  log likelihood = -5956.39
AIC=11928.79   AICc=11928.93   BIC=11968.22

Training set error measures:
                    ME     RMSE     MAE MPE MAPE      MASE       ACF1
Training set 0.1506172 82.50546 47.1043 NaN  Inf 0.7224607 0.01707355
Code
summary(midwest_arima)
Series: midwest_train_ts 
ARIMA(5,0,3) with zero mean 

Coefficients:
         ar1      ar2      ar3      ar4      ar5      ma1     ma2      ma3
      0.1952  -0.6269  -0.3389  -0.2999  -0.4563  -1.2016  0.9332  -0.3022
s.e.  0.0607   0.0312   0.0356   0.0276   0.0442   0.0702  0.0718   0.0382

sigma^2 = 21712:  log likelihood = -6544.65
AIC=13107.31   AICc=13107.49   BIC=13151.67

Training set error measures:
                    ME     RMSE      MAE MPE MAPE      MASE       ACF1
Training set 0.3936445 146.7706 85.18588 NaN  Inf 0.3810646 -0.0336952
Code
summary(south_arima)
Series: south_train_ts 
ARIMA(5,0,4) with zero mean 

Coefficients:
         ar1      ar2      ar3      ar4      ar5      ma1     ma2      ma3
      0.1881  -0.9049  -0.0800  -0.4542  -0.5567  -0.9401  1.0118  -0.7340
s.e.  0.0349   0.0292   0.0464   0.0251   0.0303   0.0357  0.0434   0.0525
         ma4
      0.4897
s.e.  0.0402

sigma^2 = 40663:  log likelihood = -6864.82
AIC=13749.64   AICc=13749.86   BIC=13798.93

Training set error measures:
                    ME     RMSE      MAE MPE MAPE      MASE        ACF1
Training set 0.1422268 200.7611 119.8731 NaN  Inf 0.3732984 -0.06523796
Code
summary(west_arima)
Series: west_train_ts 
ARIMA(4,0,3) with zero mean 

Coefficients:
          ar1     ar2      ar3      ar4      ma1      ma2     ma3
      -0.0876  0.1258  -0.5602  -0.4779  -0.5985  -0.3463  0.5343
s.e.   0.0505  0.0450   0.0330   0.0309   0.0509   0.0784  0.0421

sigma^2 = 11991:  log likelihood = -6241.45
AIC=12498.9   AICc=12499.04   BIC=12538.33

Training set error measures:
                    ME     RMSE      MAE MPE MAPE      MASE        ACF1
Training set 0.1545496 109.1261 72.50105 NaN  Inf 0.4607859 -0.08454213
Code
## Forecasting Test Data -------------------------------------------------------
east_values <- forecast(east_arima, h = length(east_test_ts))
midwest_values <- forecast(midwest_arima, h = length(midwest_test_ts))
south_values <- forecast(south_arima, h = length(south_test_ts))
west_values <- forecast(west_arima, h = length(west_test_ts))

8.2 Accuracy

Code
## Plotting & Accuracy

## East Region
east_auto_acc <- forecast::accuracy(east_values)

east_auto_acc  %>% 
  kbl(caption = "<center>Auto ARIMA East Accuracy<center>") %>% 
  kable_styling()
Auto ARIMA East Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1506172 82.50546 47.1043 NaN Inf 0.7224607 0.0170736
Code
east_auto_plot <- forecast::autoplot(east_values) +
  geom_line(data = east_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "East Region Auto ARIMA",
       y = "Daily Deaths")

east_auto_plot 

Code
## Midwest Region
midwest_auto_acc <- forecast::accuracy(midwest_values)

midwest_auto_acc %>% 
  kbl(caption = "<center>Auto ARIMA Midwest Accuracy<center>") %>% 
  kable_styling()
Auto ARIMA Midwest Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.3936445 146.7706 85.18588 NaN Inf 0.3810646 -0.0336952
Code
midwest_auto_plot <- forecast::autoplot(midwest_values) +
  geom_line(data = midwest_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "Midwest Region Auto ARIMA",
       y = "Daily Deaths")

midwest_auto_plot 

Code
## South Region
south_auto_acc <- forecast::accuracy(south_values)

south_auto_acc %>% 
  kbl(caption = "<center>Auto ARIMA South Accuracy<center>") %>% 
  kable_styling()
Auto ARIMA South Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1422268 200.7611 119.8731 NaN Inf 0.3732984 -0.065238
Code
south_auto_plot <- forecast::autoplot(south_values) +
  geom_line(data = south_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "South Region Auto ARIMA",
       y = "Daily Deaths")

south_auto_plot

Code
## West Region
west_auto_acc <- forecast::accuracy(west_values)

west_auto_acc %>% 
  kbl(caption = "<center>Auto ARIMA West Accuracy<center>") %>% 
  kable_styling()
Auto ARIMA West Accuracy
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1545496 109.1261 72.50105 NaN Inf 0.4607859 -0.0845421
Code
west_auto_plot <- forecast::autoplot(west_values) +
  geom_line(data = west_test_plot, aes(x = num, y = daily_deaths)) +
  scale_x_continuous(limits = c(950,1136),
                     breaks = c(977, 1023, 1069, 1115),
                     labels = c("Sept 2022", "Nov 2022", "Jan 2023", "March 2023")) +
  labs(title = "West Region Auto ARIMA",
       y = "Daily Deaths")

west_auto_plot

Code
# saving results
save(east_auto_plot,
     midwest_auto_plot,
     south_auto_plot,
     west_auto_plot,
     east_auto_acc,
     midwest_auto_acc,
     south_auto_acc,
     west_auto_acc,
     file = "results/uni_autoarima.rda")

9. Compare Arima Models

Code
load("results/uni_arima.rda")
load("results/uni_sarima.rda")
load("results/uni_autoarima.rda")

East Results

Code
east1 <- east_initial_arima_acc %>% 
  as.data.frame()

east2 <- east_tuned_arima_acc %>% 
  as.data.frame()

east3 <- east_initial_sarima_acc%>% 
  as.data.frame() 

east4 <- east_tuned_sarima_acc%>% 
  as.data.frame()

east5 <- east_auto_acc%>% 
  as.data.frame()

tribble(
  ~model, ~mae, ~mase,
  #--|--|----
  "initial arima", east1$MAE, east1$MASE,
  "tuned arima", east2$MAE, east2$MASE,
  "initial sarima", east3$MAE, east3$MASE,
  "tuned sarima", east4$MAE, east4$MASE,
  "auto arima", east5$MAE, east5$MASE,
) %>% 
  kbl() %>% 
  kable_styling()
model mae mase
initial arima 64.66739 0.9918341
tuned arima 47.18823 0.7237480
initial sarima 61.27670 0.9398295
tuned sarima 53.94287 0.8273471
auto arima 47.10430 0.7224607

Midwest Results

Code
midwest1 <- midwest_initial_arima_acc %>% 
  as.data.frame()

midwest2 <- midwest_tuned_arima_acc %>% 
  as.data.frame()

midwest3 <- midwest_initial_sarima_acc%>% 
  as.data.frame() 

midwest4 <- midwest_tuned_sarima_acc%>% 
  as.data.frame()

midwest5 <- midwest_auto_acc%>% 
  as.data.frame()

tribble(
  ~model, ~mae, ~mase,
  #--|--|----
  "initial arima", midwest1$MAE, midwest1$MASE,
  "tuned arima", midwest2$MAE, midwest2$MASE,
  "initial sarima", midwest3$MAE, midwest3$MASE,
  "tuned sarima", midwest4$MAE, midwest4$MASE,
  "auto arima", midwest5$MAE, midwest5$MASE,
) %>% 
  kbl() %>% 
  kable_styling()
model mae mase
initial arima 138.93519 0.6215031
tuned arima 95.00769 0.4250009
initial sarima 138.93519 0.6215031
tuned sarima 91.79687 0.4106378
auto arima 85.18588 0.3810646

South Results

Code
south1 <- south_initial_arima_acc %>% 
  as.data.frame()

south2 <- south_tuned_arima_acc %>% 
  as.data.frame()

south3 <- south_initial_sarima_acc%>% 
  as.data.frame() 

south4 <- south_tuned_sarima_acc%>% 
  as.data.frame()

south5 <- south_auto_acc%>% 
  as.data.frame()

tribble(
  ~model, ~mae, ~mase,
  #--|--|----
  "initial arima", south1$MAE, south1$MASE,
  "tuned arima", south2$MAE, south2$MASE,
  "initial sarima", south3$MAE, south3$MASE,
  "tuned sarima", south4$MAE, south4$MASE,
  "auto arima", south5$MAE, south5$MASE,
) %>% 
  kbl() %>% 
  kable_styling()
model mae mase
initial arima 203.7282 0.6344327
tuned arima 120.2240 0.3743912
initial sarima 203.7282 0.6344327
tuned sarima 120.2240 0.3743912
auto arima 119.8731 0.3732984

West Results

Code
west1 <- west_initial_arima_acc %>% 
  as.data.frame()

west2 <- west_tuned_arima_acc %>% 
  as.data.frame()

west3 <- west_initial_sarima_acc%>% 
  as.data.frame() 

west4 <- west_tuned_sarima_acc%>% 
  as.data.frame()

west5 <- west_auto_acc%>% 
  as.data.frame()

tribble(
  ~model, ~mae, ~mase,
  #--|--|----
  "initial arima", west1$MAE, west1$MASE,
  "tuned arima", west2$MAE, west2$MASE,
  "initial sarima", west3$MAE, west3$MASE,
  "tuned sarima", west4$MAE, west4$MASE,
  "auto arima", west5$MAE, west5$MASE,
) %>% 
  kbl() %>% 
  kable_styling()
model mae mase
initial arima 101.08083 0.6424269
tuned arima 63.93076 0.4063168
initial sarima 101.08083 0.6424269
tuned sarima 71.27948 0.4530221
auto arima 72.50105 0.4607859